function [psi,phi] = cmsc(psi,phi,Q,in,mu,w,dx,g)
% function [psi,phi] = cmsc(psi,phi,Q,in,mu,w,dx)
%   Inputs:
%           psi     -- fine mesh angular flux
%           phi     -- fine mesh scalar flux
%           Q       -- fine mesh external source
%           in      -- inpute structure
%           mtt     -- material index for each fine mesh
%           mu      -- angle set
%           wt      -- angle weight
%           dx      -- fine mesh delta x
% 
%   This function implements a coarse mesh, transport-based acceleration
%   scheme.  Like CMFD, we use homogenization of partially converged fine
%   mesh solution coupled with a discontinuity factor to reproduce aspects
%   of the global solution.  Following the paper of Grassi (NSE 155, 208-
%   222, 2007),  discontinuity factors are used on the outgoing and 
%   incoming partial currents to preserve net current continuity between
%   coarse meshes.  The algorithm is briefly described as follows.
%
%   1)  Perform a fine mesh sweep, yielding psi_f, the fine mesh angular
%       flux.
%   2)  Homogenize cross-sections and sources via flux-volume weighting.
%   3)  From the fine mesh sweep, we have incoming (-) and outgoing (+)
%       partial currents J- and J+ for each coarse mesh "i", which become 
%       the reference values.  Now, using the homogenized coarse mesh data,
%       perform a coarse mesh sweep, from which we obtain J*- and J*+, the
%       coarse mesh values.
%   4)  Due to homogenization, the coarse mesh J's will not match the fine
%       mesh J's, and so we define f-/+(i) = J-/+(i) / J*-/+(i).  These are
%       constant for all further coarse mesh sweeps.
%   5)  For this and all other sweeps, update the coarse mesh scalar flux
%       via neutron balance:
%             f+(i) J*+(i) - f-(i) J*-(i) + Sig(i)phi(i)V(i) = Q(i)V(i)
%       and then the scattering source
%             Q(i) <-- 0.5*S(i) + 0.5*phi(i)SigS(i)
%       (where SigS was properly flux-volume averaged).
%   6)  Perform sweeps and step 5) until coarse mesh flux is converged.  
%   7)  Update psi_f' = psi_f * psi_c' / psi_c, where the prime indicates a
%       new value
% 
%   Example discretization and use
% 
%   |              |              |              |
%   |    |    |    |    |    |    |    |    |    |   <-- fine mesh
%   |              |              |              |
%  1/2            3/2            5/2            7/2  <-- coarse mesh
%

% allocate
sigt  = zeros( length(in.xfm), 1 );   % homogenized total cross-section
sigs  = zeros( length(in.xfm), 1 );   % homogenized scatter cross-section
s_c   = zeros( length(in.xfm), 1 );   % volume-averaged external source
phi_c = zeros( length(in.xfm), 1 );   % coarse mesh scalar flux
psi_c = zeros( length(in.xcm), 2 );   % coarse mesh edge angular flux
Jp    = zeros( length(in.xcm), 1 );   % coarse mesh edge rightward (NOT J+ as in description)
Jm    = zeros( length(in.xcm), 1 );   % coarse mesh edge leftward (NOT J- as in description)
Jref  = zeros( length(in.xfm), 2 );   % reference coarse mesh total incoming and outgoing partial current
Jhom  = Jref;                         % homogenized coarse mesh total p cur
f     = Jref;                         % discontinuity factor for incoming and outgoingpartial currents
q_c   = zeros( length(in.xfm), 1 );   % coarse mesh emission density

% coarse mesh delta x's
h = in.xcm(2:end)-in.xcm(1:end-1);

s = zeros(sum(in.xfm),1);  % the angle-integrated ext src (i.e ausume isotropic)
for i = 1:length(s)
    s(i) = sum( Q(i,:,1)'.*w(:) );
end
Jp(1) = partcur( 1, psi, 1, mu, w, in );
Jm(1) = partcur(-1, psi, 1, mu, w, in );
for i = 1:length(in.xfm)  
    % note, the homogenization is meaningless until I implement a way to
    % specify coarse meshes that contain heterogeneities
    idx1     = 1 + sum( in.xfm(1:(i-1)) );            % lower index
    idx2     = sum( in.xfm(1:(i  )) );                % upper index   
    Jp(i+1)  = partcur( 1, psi, idx2+1, mu, w, in );  % J+ -->
    Jm(i+1)  = partcur(-1, psi, idx2+1, mu, w, in );  % J- <--
    Jref(i,1)= Jp(i)+Jm(i+1);
    Jref(i,2)= Jm(i)+Jp(i+1);
    phi_c(i) = sum( dx(idx1:idx2)'.*phi(idx1:idx2,1) ) / h(i); % coarse mesh phi
    sigt(i)  = sum( dx(idx1:idx2)'.*phi(idx1:idx2,1)*in.data( in.mt(i), 1 ));
    sigt(i)  = sigt(i)/(phi_c(i)*h(i)); % total cross-section
    sigs(i)  = sum( dx(idx1:idx2)'.*phi(idx1:idx2,1)*in.data( in.mt(i), 5 ));    
    sigs(i)  = sigs(i)/(phi_c(i)*h(i));
    s_c(i)   = sum( dx(idx1:idx2)'.*s(idx1:idx2,1) ) / h(i);    
end

phi_ref = phi_c;

% compute Qhat

% here, we use an S2-like quadrature.  We want to conserve partial
% currents.  In S2, we have just one angle going right and one going left.
% The partial current J+ = w*mu*psi, or psi = J+/w/mu, and likewise for the
% leftward J-.  We correct Q such that psi(i+1)=J+/w/mu=f(psi(i),Q).

% using an S2-like quadrature
mu_c = 0.01;%0.5773502691;
wt_c = 0.5;

% % going right
% psi_c(1,1) = 0; % vacuum
% for i = 1:length(in.xfm)
%     psi_c(i+1,1) = Jp(i+1)/mu_c/wt_c;            % reference right edge psi
%     ef           = exp( -h(i)*sigt(i)/mu_c );    % exponential factor
%     % want psi2 = psi2ref
%     %  psi2ref = 0.5*(Qsc + Qex + Qhat)/sigt * (1-ef) + psi1*ef
%     %  Qhat = (psi2ref-psi1*ef) * sigt/0.5 - Qsc -Qex
%     Qhat(i,1) = (psi_c(i+1,1)-psi_c(i,1)*ef)*sigt(i)/(1-ef)*2 - phi_c(i)*sigs(i)-s_c(i);
%     %tmpPsi2 = 0.5*(phi_c(i)*sigs(i)+s_c(i)+Qhat(i))/sigt(i) * (1-ef) + psi_c(i,1)*ef;
% end
% % going left
% psi_c(end,2) = 0;
% for i = length(in.xfm):-1:1
%     psi_c(i,2) = Jm(i)/mu_c/wt_c;            % reference right edge psi
%     ef         = exp( -h(i)*sigt(i)/mu_c );    % exponential factor
%     Qhat(i,2)  = (psi_c(i,2) - psi_c(i+1,2)*ef)*sigt(i)/(1-ef)*2 - phi_c(i)*sigs(i)-s_c(i);
%     % check
%     %tmpPsi2 = 0.5*(phi_c(i)*sigs(i)+s_c(i)+Qhat(i))/sigt(i) * (1-ef) + psi1*ef;
% end
% 
% Qhat
% Qhat(:,:) = 0;
% hold on
% plot(x,phi,xc,phi_c,'o')
% poop = 1;
% phi_c
% phi_c = scalarflux( psi_c, s_c, phi_c.*sigs, Qhat, sigt, h )

% We have reference partial currents going left and right
% convergence parameters
eps_phi = 1e-6; 
max_it  = 1000;
err_phi = 1;    
it = 0;
% Begin coarse mesh source iterations
NCM = length(phi_c);

ef = zeros(NCM,1);
for i = 1:NCM       % left-to-right
    ef(i)  = exp(-h(i)*sigt(i)/mu_c);
end
while (err_phi > eps_phi && it <= max_it )
    % Save old scalar flux
    phi0 = phi_c; 
    % Update sources
    for i = 1:NCM
        q_c(i) = s_c(i) + sigs(i)*phi_c(i);
    end
    % Perform sweeps
    for i = 1:NCM       % left-to-right
        psi_c(i+1,1) = 0.5*(q_c(i))/sigt(i)*(1-ef(i)) + psi_c(i,1)*ef(i);
    end
    for i = NCM:-1:1  	% left-to-right
        psi_c(i,2)   = 0.5*(q_c(i))/sigt(i)*(1-ef(i)) + psi_c(i+1,2)*ef(i); 
    end    
    for i = 1:NCM
        Jp(i+1)      = psi_c(i+1,1) * mu_c * wt_c;
        Jm(NCM+1-i)  = psi_c(NCM+1-i,2) * mu_c * wt_c;
    end
    for i = 1:NCM
        Jhom(i,1)= Jp(i)+Jm(i+1);
        Jhom(i,2)= Jm(i)+Jp(i+1);        
    end
    if ( it == 0 )
       f = Jref./Jhom;
    end
    % Update phi from neutron balance
    for i = 1:NCM
       phi_c(i) = ( q_c(i)*h(i) + f(i,1)*Jhom(i,1) - f(i,2)*Jhom(i,2) ) ...
                  / ( sigt(i)*h(i) );
    end
    % Update error and iteration counter
    err_phi =  max(  abs(phi_c-phi0)./phi0 );
    it = it + 1;
end

for i = 1:length(in.xfm)  
    % within cell indices
    idx1 = 1 + sum( in.xfm(1:(i-1)) ); % lower index
    idx2 = sum( in.xfm(1:(i  )) );     % upper index   
    phi(idx1:idx2,1)   = phi(idx1:idx2,1)  * phi_c(i)/phi_ref(i);
end
%figure(2),plot(xc,phi_ref,'k',xc,phi_c,'r--')
%lala=1;
end

function J = partcur( flag, psi, i, mu, w, in )
    if ( flag == 1 )
        idxs = in.ord/2+1:in.ord;
    else
        idxs = 1:in.ord/2;
    end
    J = sum( abs(mu(idxs)) .* w(idxs) .* psi(i,idxs)' );
end

